///@file bergamote.cpp
///@brief Bergamote main program
///@author Arnaud Duval
///@version 0.0
///@date 2010/05/10

///@mainpage My Personal Index Page
///@section intro_sec Introduction
///This is the introduction.
///@section install_sec Installation
///@subsection step1 Step 1: Opening the box
///etc...



#include <iostream>
#include <fstream>
#include <cstdlib>
#include "config.h"
#include "mathobject.h"
#include "FEProblem.h"
#include "random/include/IndexSet.h"
#include "random/include/ProbUtils.h"
#include "random/include/SobolSequence.h"
#include "signal/include/Sample.h"
#include "StrainTensor.h"
#include "StressTensor.h"
#include "results/include/ResultTable.h"
#include <ctime>
//#include "gui/include/MainWindow.h"
#include "tests/include/PerfTests.h"
#include "utils/include/StringUtils.h"
#include "utils/include/FileUtils.h"
#include "tests/include/MathobjectTests.h"
#include "ElastoplasticMaterial.h"
#include <cmath>

///@brief print title art
void TitleArt();
///@brief Parse Command line
void CommandParser(int, char**);


///@brief Test routine for polynomial chaos expansion
void TestCP();
///@brief Test routine for UMAT in 3D
void TestUMAT3D();
///@brief Test routine for UMAT in 2D
void TestUMAT2D();
///@brief Function to debug a 2D UMAT
void DebugUMAT2D();



int main(int argc, char **argv)
{


    
    //DebugUMAT2D();
    TestUMAT2D();
  
    //TitleArt();
    //TestVector();
    /*Mesh monMaillage;
    monMaillage.ReadXML("/home/arnaud/plaque.xml");
    //std::cout << monMaillage << "\n";
    Quad4 test;
    FEProblem houston(&monMaillage);
  
    houston.AddNodeGroup("blocage");
    houston.AddNodeGroup("Déplacement imposé");
    houston.AddNodeToGroup("n_1", "blocage");
    houston.AddNodeToGroup("n_2", "blocage");
    houston.AddNodeToGroup("n_3", "blocage");
    houston.AddNodeToGroup("n_7", "Déplacement imposé");
    houston.AddNodeToGroup("n_8", "Déplacement imposé");
    houston.AddNodeToGroup("n_9", "Déplacement imposé");*/


	///Ci-dessous : voir à quoi ça sert et s'il faut le garder ou pas
	/*
	unsigned int n = 100000;	//Nombre de réalisations
	Distribution F(n, "uniform", 48., 52.);
	Distribution L(n, "gauss", 1000., 10.);
	Distribution b(n, "gauss", 20., 1.);
	Distribution h(n, "gauss", 30., 1.5);
	Distribution age(n, "uniform", 55., 75.);
	Distribution E(n, "gauss", 210000., 10000.);
	//Ici, les distributions sont générées.
	std::ofstream f("f.txt");
	f << "none" << "\t" << 0. << "\t" << 0. << "\t" << n << "\n";
	for(unsigned int i = 1 ; i <= n ; i++)
	{
		double rF = F[i];
		double rL = L[i];
		double rb = b[i];
		double rh = h[i];
		double rE = E[i];
		double res = 4.*rF*pow(rL, 3.)/(rE*rb*pow(rh, 3.));
		f << res << "\n";
	}
	f.close();

	Distribution fleche;
	fleche.ReadFromFile("f.txt");
	std::cout << "Moyenne : " << fleche.Average() << "\n";
	std::cout << "Ecart-type : " << fleche.StandardDeviation() << "\n";
	double *c95;
	double *c99;
	c95 = fleche.Confidence(95.);
	c99 = fleche.Confidence(99.);
	std::cout << "Intervalle de confiance 95 % : " << c95[0] << "\t" << c95[1] << "\n";
	std::cout << "Intervalle de confiance 99 % : " << c99[0] << "\t" << c99[1] << "\n";

	Histogram h1(fleche, 10);
	Histogram h2(fleche, 100);
	Histogram h3(fleche, 1000);
	Histogram h4(fleche, 10000);

	std::ofstream histo10("histo10.txt");
	std::ofstream histo100("histo100.txt");
	std::ofstream histo1000("histo1000.txt");
	std::ofstream histo10000("histo10000.txt");

	histo10 << h1;
	histo100 << h2;
	histo1000 << h3;
	histo10000 << h4;

	histo10.close();
	histo100.close();
	histo1000.close();
 	histo10000.close();

	return 0;
	
	*/




	//CommandParser(argc, argv);
	return 0;

	//Test of PCExpansion class
	/*PCExpansion pce("/home/arnaud/Desktop/poutre/PCExpansion.pce");
	std::cout << pce << "\n";
	std::ifstream inF("/home/arnaud/Desktop/poutre/F.txt");
	std::ifstream inL("/home/arnaud/Desktop/poutre/L.txt");
	std::ifstream inb("/home/arnaud/Desktop/poutre/b.txt");
	std::ifstream inh("/home/arnaud/Desktop/poutre/h.txt");
	std::ifstream inage("/home/arnaud/Desktop/poutre/age.txt");
	std::ifstream inE("/home/arnaud/Desktop/poutre/E.txt");
	std::ofstream res("/home/arnaud/Desktop/poutre/res-pce.txt");

	double F,L,b,h,A,E;

	double reals[100000];

	for(int i = 0 ; i < 100000 ; i++)
	{
		inF >> F;
		inL >> L;
		inb >> b;
		inh >> h;
		inage >> A;
		inE >> E;
		Vector real(6);
		real[1] = F;
		real[2] = L;
		real[3] = b;
		real[4] = h;
		real[5] = A;
		real[6] = E;
		res.precision(18);
		reals[i] = pce.Evaluate(real);
		res << reals[i] << "\n";
		//stat[(int)(f*100.)] ++;
	}

	res.close();
	inF.close();
	inL.close();
	inb.close();
	inh.close();
	inage.close();
	inE.close();


	Histogram stat(reals, 100000, 100);

	std::ofstream statres("/home/arnaud/Desktop/poutre/stat-pce.txt");
	statres << stat;

	statres.close();*/




	/*Distribution F,L,b,h,age,E;
	F.ReadFromFile("distributions/F.txt");
	L.ReadFromFile("distributions/L.txt");
	b.ReadFromFile("distributions/b.txt");
	h.ReadFromFile("distributions/h.txt");
	age.ReadFromFile("distributions/age.txt");
	E.ReadFromFile("distributions/E.txt");*/

	
	//std::cout << d ;
	//Histogram h(d, 100);		//Erreur de seg ici
	//std::cout << h << "\n";
	//std::cout << d.Average() << "\n";
	//std::cout << d.StandardDeviation() << "\n";
	//d.BubbleSort();
	//std::cout << d << "\n";
	/*for(unsigned int i = 10 ; i < 100000 ; i*=2)
	{
		Distribution d(i, "gauss", 10000., 1000.);
		double *a=d.Confidence(95.);
		std::cout << i  << "\t" << a[0] << "\t" << a[1] << "\n";
	}
	return 0;*/


	//IndexSet res = test.I_i(2);
	//std::cout << "#######\n";
	//std::cout << res << "\n";



	return 0;

	//******Moments quadratiques, cas-test du tube de section rectangulaire (60*30 ext 55*25 int)

	/*Point a(2);
	Point b(2);
	Point c(2);
	Point d(2);
	Point e(2);
	Point f(2);
	Point g(2);
	Point h(2);
	Point i(2);
	Point j(2);
	Point k(2);
	Point l(2);

	a.Coord(1) = 15.;
	a.Coord(2) = 30.;
	b.Coord(1) = 15.;
	b.Coord(2) = -30.;
	c.Coord(1) = -15.;
	c.Coord(2) = -30.;
	d.Coord(1) = -15.;
	d.Coord(2) = 30.;
	e.Coord(1) = -12.5;
	e.Coord(2) = 27.5;
	f.Coord(1) = 12.5;
	f.Coord(2) = 27.5;
	g.Coord(1) = 12.5;
	g.Coord(2) = -27.5;
	h.Coord(1) = -12.5;
	h.Coord(2) = -27.5;
	i.Coord(1) = 0.;
	i.Coord(2) = 30.;
	j.Coord(1) = 0.;
	j.Coord(2) = -30.;
	k.Coord(1) = 15.;
	k.Coord(2) = 0.;
	l.Coord(1) = -15.;
	l.Coord(2) = 0.;

	Triangle t1(e,i,d);
	Triangle t2(e,f,i);
	Triangle t3(f,a,i);
	Triangle t4(f,k,a);
	Triangle t5(g,k,f);
	Triangle t6(g,b,k);
	Triangle t7(j,b,g);
	Triangle t8(j,g,h);
	Triangle t9(j,h,c);
	Triangle t10(c,h,l);
	Triangle t11(h,e,l);
	Triangle t12(l,e,d);

	double sumxx = 0.;
	double sumyy = 0.;
	sumxx += t1.Ixx_O();
	sumxx += t2.Ixx_O();
	sumxx += t3.Ixx_O();
	sumxx += t4.Ixx_O();
	sumxx += t5.Ixx_O();
	sumxx += t6.Ixx_O();
	sumxx += t7.Ixx_O();
	sumxx += t8.Ixx_O();
	sumxx += t9.Ixx_O();
	sumxx += t10.Ixx_O();
	sumxx += t11.Ixx_O();
	sumxx += t12.Ixx_O();

	sumyy += t1.Iyy_O();
	sumyy += t2.Iyy_O();
	sumyy += t3.Iyy_O();
	sumyy += t4.Iyy_O();
	sumyy += t5.Iyy_O();
	sumyy += t6.Iyy_O();
	sumyy += t7.Iyy_O();
	sumyy += t8.Iyy_O();
	sumyy += t9.Iyy_O();
	sumyy += t10.Iyy_O();
	sumyy += t11.Iyy_O();
	sumyy += t12.Iyy_O();

	std::cout << sumxx << "\n";
	std::cout << sumyy << "\n";*/



	/*ResultTable t("/home/arnaud/Desktop/resultats-2plaques/real_0008_depX_163.freq", 64, 2, false);
	double ** maxs;
	unsigned int count;
	maxs = t.GetMaxLines(2, count);
	for(int i = 0 ; i < count ; i++)
		std::cout << maxs[0][i] << "\t" << maxs[1][i] << "\n";
	

	for(int i = 0 ; i < 2 ; i++)
		delete[] maxs[i];
	delete[] maxs;*/
	
	

	/*std::string str[4];
	str[0] = "gauss";
	str[1] = "uniform";
	str[2] = "gauss";
	str[3] = "gauss";

	PCbasis test(4, str, 5);
	std::cout << test << "\n";*/



	//TestMatrix(10000);

	//TestCP();
	return 0;

	/*std::cout << FileExists("/home/arnaud/result_P1-1pc.res") << "\n";
	std::cout << FileExists("/home/arnaud/result_P1-2pc.res") << "\n";
	std::cout << FileExists("/home/arnaud/result_P1-3pc.res") << "\n";

	ResultTable table1("/home/arnaud/result_P1-1pc.res", 40969u, 5u, false);
	ResultTable table2("/home/arnaud/result_P1-2pc.res", 40969u, 5u, false);
	ResultTable table3("/home/arnaud/result_P1-3pc.res", 40969u, 5u, false);

	double *maxline;
	maxline = table1.GetMinLine(2u);
	std::cout << "1 pourcent : " << maxline[0] << "\t" << maxline[1] << "\n";
	maxline = table2.GetMinLine(2u);
	std::cout << "2 pourcent : " << maxline[0] << "\t" << maxline[1] << "\n";
	maxline = table3.GetMinLine(2u);
	std::cout << "3 pourcent : " << maxline[0] << "\t" << maxline[1] << "\n";
	std::cout << "\n\n";
	maxline = table1.GetMinLine(3u);
	std::cout << "1 pourcent : " << maxline[0] << "\t" << maxline[2] << "\n";
	maxline = table2.GetMinLine(3u);
	std::cout << "2 pourcent : " << maxline[0] << "\t" << maxline[2] << "\n";
	maxline = table3.GetMinLine(3u);
	std::cout << "3 pourcent : " << maxline[0] << "\t" << maxline[2] << "\n";
	std::cout << "\n\n";
	maxline = table1.GetMinLine(4u);
	std::cout << "1 pourcent : " << maxline[0] << "\t" << maxline[3] << "\n";
	maxline = table2.GetMinLine(4u);
	std::cout << "2 pourcent : " << maxline[0] << "\t" << maxline[3] << "\n";
	maxline = table3.GetMinLine(4u);
	std::cout << "3 pourcent : " << maxline[0] << "\t" << maxline[3] << "\n";


	//return 0;

	//TestMatrix(1000);

	

	//**********Test lecture maillage
	//Mesh monmaillage;
	//monmaillage.ReadXML("essai.xml");
	//monmaillage.Equilibrate();
	//std::cout << monmaillage << "\n";

	//**********Test *********************
	StressTensor s(3);
	s(1,1) = 240.;
	s(1,2) = 127.;
	s(1,3) = 218.;
	s(2,3) = 47.;
	s(2,2) = 25.;
	s(3,3) = 21.;

	StrainTensor e;
	e &= s;
	//std::cout << e << "\n";
	//std::cout << s << "\n";
	//std::cout << (Tensor2)e << "\n";
	

	

	//**********Test stress tensor*********
	//StressTensor t(3);
	//t(1,1) = 240.;
	//t(1,2) = 127.;
	//t(1,3) = 218.;
	//t(2,3) = 47.;
	//t(2,2) = 25.;
	//t(3,3) = 21.;
	//std::cout << t << "\n";
	//std::cout << t.EquivVM() << "\n";

	//**********Test strain tensor*********
	/*StrainTensor t(3);
	t(1,1) = 0.05;
	t(2,2) = -0.02;
	t(3,3) = -0.03;
	t(1,2) = 0.045;
	t(1,3) = -0.01;
	t(2,3) = 0.02;

	std::cout << t << "\n";
	std::cout << t.EquivVM() << "\n";
	std::cout << (StrainVector)t << "\n";

	StrainTensor a(2);
	a(1,1) = 0.03;
	a(2,2) = -0.01;
	a(1,2) = 0.002;

	std::cout << a << "\n";
	std::cout << a.EquivVM() << "\n";
	std::cout << (StrainVector)a << "\n";*/
	


	
	//***********Test polynômes************
	/*Polynomial p(6, 'H');
	std::cout << p.Evaluate(23.6) << "\n";*/

	//SobolSequence a(3, 25);
	//std::cout << a << "\n";



	
	
	


	return 0;
}
/*
///@brief print title art
void TitleArt()
{
	std::cout << " _                                          _        \n";
	std::cout << "| |                                        | |       \n";
	std::cout << "| |__   ___ _ __ __ _  __ _ _ __ ___   ___ | |_ ___  \n";
	std::cout << "| '_ \\ / _ \\ '__/ _` |/ _` | '_ ` _ \\ / _ \\| __/ _ \\ \n";
	std::cout << "| |_) |  __/ | | (_| | (_| | | | | | | (_) | ||  __/ \n";
	std::cout << "|_.__/ \\___|_|  \\__, |\\__,_|_| |_| |_|\\___/ \\__\\___| \n";
	std::cout << "                 __/ |                               \n";
	std::cout << "                |___/                                \n";
	std::cout << "bergamote " << BERGAMOTE_MAJOR<< "." << BERGAMOTE_MINOR << "." << BERGAMOTE_PATCH << "\n";
}

///@brief Parse Command line
///@param n number of words in the command line
///@param a array of words of the command line
void CommandParser(int n, char** a)
{
	if(n > 1)
	{
		std::string action = a[1];
		if(action == "-fft")
		{
			if(n != 4)
			{
				if(n == 2)
				{
					std::cerr << "Error : missing input and output file. Usage : bergamote -fft infile outfile\n";
					std::cerr << "Exiting...\n";
				}
				if(n == 3)
				{
					std::cerr << "Error : missing output file. Usage : bergamote -fft infile outfile\n";
					std::cerr << "Exiting...\n";
				}
				return;
			}
			if(!FileExists(a[2]))
			{
				std::cerr << "Error : unable to open file " << a[2] << "\n";
				std::cerr << "Exiting...\n";
				return;
			}
			Sample s;
			s.Read(a[2]);
			s.WriteFFT(a[3]);
			std::cout << "Output written to " << a[3] << "\n";
		}
		else if(action == "-post")
		{
			if(n != 3)
			{
				std::cerr << "Error : missing input file. Usage : bergamote -post infile\n";
				std::cerr << "Exiting...\n";
				return;
			}
			if(!FileExists(a[2]))
			{
				std::cerr << "Error : unable to open file " << a[2] << "\n";
				std::cerr << "Exiting...\n";
				return;
			}
			unsigned int nblines, nbcols;
			std::cout << "\nNumber of columns to read :";
			std::cin >> nbcols;
			std::cout << "\nNumber of lines to read :";
			std::cin >> nblines;
			ResultTable table(a[2], nblines, nbcols, false);
			std::cout << "\nNumber of columns in the output file :";
			unsigned int nbcolsout;
			unsigned int *indexout;
			std::cin >> nbcolsout;
			indexout = new unsigned int[nbcolsout];
			for(int i = 0 ; i < nbcolsout ; i++)
			{
				std::cout << "\nindex of columns in the input file to be columns nb " << i+1 << " in the out file :";
				std::cin >> indexout[i];
			}
			std::cout << "\nName of the output file :";
			std::string outfilename;
			std::cin >> outfilename;
			table.WriteResult(outfilename, nbcolsout ,indexout);
			std::cout << "Output written to " << outfilename << "\n";
		}
		else if(action == "-chaos")
		{
			if(n != 3)
			{
				std::cerr << "Error : missing input file. Usage : bergamote -chaos definitionfile.chaos\n";
				std::cerr << "Exiting...\n";
				return;
			}
			if(!FileExists(a[2]))
			{
				std::cerr << "Error : unable to open file " << a[2] << "\n";
				std::cerr << "Exiting...\n";
				return;
			}
			ChaosPrepro(a[2]);
		}
		else if(action == "-chaos-sol")
		{
			if(n != 3)
			{
				std::cerr << "Error : missing input file. Usage : bergamote -chaos-sol definitionfile.chaos-sol\n";
				std::cerr << "Exiting...\n";
				return;
			}
			if(!FileExists(a[2]))
			{
				std::cerr << "Error : unable to open file " << a[2] << "\n";
				std::cerr << "Exiting...\n";
				return;
			}
			ChaosSolver(a[2]);
		}
		else if(action == "-section")
		{
			if(n != 3)
			{
				std::cerr << "Error : missing input file. Usage : bergamote -section definitionfile.section\n";
				std::cerr << "Exiting...\n";
				return;
			}
			if(!FileExists(a[2]))
			{
				std::cerr << "Error : unable to open file " << a[2] << "\n";
				std::cerr << "Exiting...\n";
				return;
			}
			ClosedPolyline line(a[2]);
			CrossSection sec(line);
			std::cout << "A :\t" << sec.GetA() << "\n";
			std::cout << "Ixx :\t" << sec.GetIxx() << "\n";
			std::cout << "Iyy :\t" << sec.GetIyy() << "\n";
			std::cout << "Ixy :\t" << sec.GetIxy() << "\n";
			std::cout << "CGx :\t" << sec.GetCGx() << "\n";
			std::cout << "CGy :\t" << sec.GetCGy() << "\n";
			std::cout << "Iw :\t" << sec.GetIw() << "\n";
			std::cout << "J :\t " << sec.GetJ() << "\n";
			std::cout << "SHx :\t" << sec.GetSHx() << "\n";
			std::cout << "SHy :\t" << sec.GetSHy() << "\n";
		}
		else
		{
			return;
		}
	}
}





///@brief Test routine for polynomial chaos expansion
void TestCP()
{
	SobolSequence seq(4,50);
	//std::cout << seq << "\n";
	
	seq.TransformGauss(1,1.831e-5,0.05e-5);
	seq.TransformGauss(2,2.9e-3,0.2e-3);
	seq.TransformUniform(3,0.08,0.12);
	seq.TransformGauss(4,200.e9, 10.e9);
	std::cout << seq << "\n";
}
*/

///@brief Test routine for UMAT
void TestUMAT3D()
{

    /////TEST D UNE LOI DE COMPORTEMENT UMAT
    
    ElastoplasticMaterial mat(210000., 0.3, 450., 2., 300.);

    Vector stress(6);
    Vector strain(6);
    double statev[7];
    Vector dstrain(6);
    stress.Nullify();
    strain.Nullify();
    dstrain.Nullify();

    for(unsigned int i = 0 ; i < 7 ; i++)
    {
        statev[i] = 0.;
    }

    Matrix ddsdde(6u,6u);
    
    //Ci-dessous : version avec pilotage en contrainte

    
    double sig[1000];
    for(unsigned int idx = 0 ; idx < 1000 ; idx++)
    {
        sig[idx] = 500.*sin((double)idx / 100.);

    }

    //for(double sig = 0. ; sig <= 450. ; sig += 40.)
    for(unsigned int idx = 0 ; idx < 1000 ; idx++)
    {
        double statevini[7];
        for(unsigned int i = 0 ; i < 7 ; i++)
            statevini[i] = statev[i];
        Vector sigma_pil(6);
        sigma_pil[1] = sig[idx];
        sigma_pil[2] = 0.;
        sigma_pil[3] = 0.;
        sigma_pil[4] = 0.;
        sigma_pil[5] = 0.;
        sigma_pil[6] = 0.;
        Vector stressini = stress;
        //std::cout << "sigma_pil : "<< sigma_pil << "\n";

        Vector dstrain(6);
        dstrain.Nullify();

        //Première prédiction

        //Calcul du module tangent
        mat.UMAT(strain, dstrain, stress, ddsdde, statev, 7);

        Matrix invjac = ddsdde.Invert();
        dstrain = -1.*invjac|(stressini - sigma_pil);

        bool conv = false;

        for(unsigned int counter = 0 ; (counter <= 10)&&(!conv) ; counter++)
        {
            mat.UMAT(strain, dstrain, stress, ddsdde, statev, 7);

            if((stress - sigma_pil).Euclide_Norm() < 1.e-3)
            {
                conv = true;
            }
            else
            {
                invjac = ddsdde.Invert();
                dstrain -= (invjac|(stress - sigma_pil));
                stress = stressini;
                for(unsigned int i = 0 ; i < 7 ; i++)
                {
                    statev[i] = statevini[i];
                    stress = stressini;
                }
            }
        }
        if(!conv)
            std::cerr << "Erreur de convergence dans la programme principal\n";
        strain += dstrain;
        std::cout << strain << "\t" << stress << "\n";
    }

    
    //Ci dessous : version avec pilotage en déplacement
    /*for(double eps = 0. ; eps <= 0.1 ; eps += 0.001)
    {
        dstrain[1] = 0.001;
        dstrain[2] = -0.5*0.001;
        dstrain[3] = -0.5*0.001;
        dstrain[4] = 0.;
        dstrain[5] = 0.;
        dstrain[6] = 0.;

        mat.UMAT(strain, dstrain, stress, ddsdde, statev, 7);
        strain += dstrain;

        std::cout << strain << "\t" << stress << "\n";
    }*/


}

///@brief Test routine for UMAT
void TestUMAT2D()
{

    /////TEST D UNE LOI DE COMPORTEMENT UMAT en 2D contraintes planes

    ElastoplasticMaterial mat(210000., 0.3, 450., 2., 300.);

    Vector stress(3);
    Vector strain(3);
    double statev[4];
    Vector dstrain(3);
    stress.Nullify();
    strain.Nullify();
    dstrain.Nullify();


    for(unsigned int i = 0 ; i < 4 ; i++)
    {
        statev[i] = 0.;
    }

    Matrix ddsdde(3u,3u);

    //Ci-dessous : version avec pilotage en contrainte


    double sig[1000];
    for(unsigned int idx = 0 ; idx < 1000 ; idx++)
    {
        sig[idx] = 500.*sin((double)idx / 10.);

    }

    //for(double sig = 0. ; sig <= 450. ; sig += 40.)
    for(unsigned int idx = 0 ; idx < 1000 ; idx++)
    {
        //std::cout << "****** INCREMENT : " << idx << " *******\n";
        double statevini[4];
        for(unsigned int i = 0 ; i < 4 ; i++)
            statevini[i] = statev[i];
        Vector sigma_pil(3);
        sigma_pil[1] = 0.;//sig[idx];
        sigma_pil[2] = 0.;
        sigma_pil[3] = sig[idx];
        Vector stressini = stress;
        //std::cout << "sigma_pil : "<< sigma_pil << "\n";

        Vector dstrain(3);
        dstrain.Nullify();

        //Première prédiction

        //Calcul du module tangent
        mat.UMAT(strain, dstrain, stress, ddsdde, statev, 4);

        //std::cout << "REPERE PREDICTION\n";

        Matrix invjac = ddsdde.Invert();
        dstrain = -1.*invjac|(stressini - sigma_pil);

        bool conv = false;

        for(unsigned int counter = 0 ; (counter <= 100)&&(!conv) ; counter++)
        {
            mat.UMAT(strain, dstrain, stress, ddsdde, statev, 4);
            //std::cout << "APRES UMAT : " << stress[3] << "\t" << strain[3] + dstrain[3] << "\t" << sigma_pil[3] << "\n";
            //std::cout << "jac : \n" << ddsdde << "\n\n";

            //std::cout << "###[" << counter << "] Erreur prog principal : " << (stress - sigma_pil).Euclide_Norm() << "\n";
            //std::cout << "##Info sigma_pil[3] : " << sigma_pil[3] << "\n";
            
            if((stress - sigma_pil).Euclide_Norm() < 1.e-3)
            {
                conv = true;
            }
            else
            {
                invjac = ddsdde.Invert();
                //std::cout << "jacinv :\n" << ddsdde.Invert() << "\n\n";
                //std::cout << "corr : " << (invjac|(stress - sigma_pil)) << "\n";
                dstrain -= (invjac|(stress - sigma_pil));
                stress = stressini;
                for(unsigned int i = 0 ; i < 4 ; i++)
                {
                    statev[i] = statevini[i];
                }
            }
            //std::cout << "##Apres correct : " << "strain : " << strain << "   dstrain : " << dstrain << "\n";
        }
        if(!conv)
            std::cerr << "#######Erreur de convergence dans le programme principal\n";
        strain += dstrain;
        std::cout << strain << "\t" << stress << "\n";
    }


    //Ci dessous : version avec pilotage en déplacement
    /*double deps = 0.001;
    for(double eps = 0. ; eps <= 0.1 ; eps += deps)
    {
        dstrain[1] = 0.;//deps;
        dstrain[2] = 0.;//-0.5*deps;
        dstrain[3] = deps;

        mat.UMAT(strain, dstrain, stress, ddsdde, statev, 4);
        strain += dstrain;

        std::cout << strain << "\t" << stress << "\n";
    }*/
}

///@brief Function to debug a 2D UMAT
void DebugUMAT2D()
{
    ElastoplasticMaterial mat(210000., 0.3, 450., 2., 300.);
    Vector strain(3), dstrain(3), stress(3);
    Matrix ddsdde(3u, 3u);
    strain.Nullify();
    dstrain.Nullify();
    stress.Nullify();
    strain[3] = 4.03441;
    dstrain[3] = -0.0593845;
    stress[3] = 499.787;
    double statev[4];
    statev[0] = 2.84838;
    statev[1] = 0.;
    statev[2] = 0.;
    statev[3] = 4.02822;

    mat.UMAT(strain, dstrain, stress, ddsdde, statev, 4);
}
